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Molecular  Dynamics  of  Shocks  in  Crystals  with  Defects 


1  INTRODUCTION 

This  paper  investigates,  on  the  molecular  scale,  the  propagation  of  vibrational  distur¬ 
bances  through  various  two-dimensional  lattices,  using  the  techniques  of  computational 
non-equilibrium  molecular  dynamics.  These  disturbances  will  be  referred  to  as  “shocks”, 
as  is  common  in  the  literature,  because,  although  they  have  a  finite  width  and  a  definite 
structure,  they  are  thought  to  be  the  microscopic  analogues  of  the  travelling  discontinuities 
that  are  defined  as  shocks  in  the  macroscopic  domain  of  continuum  mechanics. 

The  motivation  behind  the  work  reported  in  this  paper  is  the  desire  to  understand 
some  of  the  peculiar  properties  of  the  initiation  and  propagation  of  detonations  in  solid 
chemical  explosives.  A  detonation,  as  the  term  is  usually  understood  in  these  contexts, 
is  a  physico-chemical  process  with  a  characteristic  structure.1  This  structure  consists, 
in  a  solid,  liquid,  or  gas,  of  a  shock  wave  or  some  other  localized,  traveling  disturbance, 
followed  by  an  associated  reaction  front,  which  separates  material  that  has  participated  in  a 
chemical  reaction  from  material  that  is  not  yet  reacted.  The  association  between  these  two 
traveling  interfaces,  mediated  by  the  “induction  zone”  of  (usually)  compressed,  stressed, 
unreacted  material  between  them,  is  a  complex  interdependence  whose  nature  has  been  the 
subject  of  many  investigations,  a  few  of  which  we  shall  discuss  below.  For  now  we  should 
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merely  point  out  that,  on  the  one  hand,  the  shock  is  the  cause  of  the  reaction  front,  for  in 
its  passage  it  creates  the  conditions  which  lead  to  a  particular  series  of  reactions,  while, 
on  the  other  hand,  the  exothermic  reaction  front  maintains  and  accelerates  the  shock  by 
supplying  kinetic  and  thermal  energy. 

In  the  fluid  phases,  those  properties  governing  the  nature  of  the  detonation  process  are 
the  chemical  composition  of  the  material  along  with  its  intensive  thermodynamic  variables, 
and,  sometimes,  the  details  of  the  method  used  to  initiate  the  detonation.  In  a  solid,  even 
when  all  of  these  quantities  are  held  constant,  there  may  still  be  significant  variability 
in  the  responses  of  a  series  of  samples.  For  example,  in  an  experiment  in  which  crystals 
are  detonated  by  dropping  a  weight  on  them,  the  variability  may  take  the  form  of  a  wide 
range  of  initial  heights  necessary  for  detonation  in  apparently  identical  samples.  A  possible 
source  of  this  variable  response  might  be  some  unknown  and  uncontrolled  variation  in  the 
experimental  conditions,  presenting  different  samples  with  significantly  different  initial 
conditions,  combined  with  an  extreme  sensitivity  to  initial  conditions  in  the  system.  But 
it  must  be  recognized  that  a  sample  of  solid  material,  unlike  a  fluid,  is  characterized  not 
only  by  its  chemical  composition  and  its  thermodynamic  state.  The  particular  spatial 
arrangement  of  molecules  constituting  each  sample  makes  them  all  unique.  The  fact  that 
the  shock,  solitary  compressional  wave,  or  other  disturbance  taking  part  in  the  structure 
of  the  detonation  process  is  a  mechanical  mode  of  vibration  supported  by  the  molecular 
lattice,  whose  nature  is  completely  dependent  on  the  detailed  structure  of  that  lattice, 
implies  that  the  spatial  arrangement  of  molecules  must  influence  the  detonation  process. 
This  much  is  uncontroversial.  What  may  be  more  surprising  is  the  suggestion2  investigated 
here,  that  certain  minute  changes  in  the  molecular  arrangement,  involving  the  locations 
of  just  a  few  molecules,  may  have  a  large,  even  determining,  effect  on  a  macroscopic 
detonation.  This  suggestion  gains  plausibility  if  one  considers  that  the  actual  thickness  of 
the  shock  front  may  be  only  a  few  lattice  planes.3  Therefore  the  transfer  of  kinetic  energy 
from  a  molecule  to  its  neighbors  in  the  direction  of  propagation  involves  the  participation 
of  only  a  few  molecules  (along  the  wavevector  direction)  in  a  particular  vibrational  mode 
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at  any  one  time. 

One  traditional  theoretical  approach  to  understanding  detonations,  and  the  initiation 
of  a  detom  tion  by  a  shock  wave,  involves  a  description  in  terms  of  thermodynamic  variables 
that  are  related  by  an  equation  of  state  for  the  material.1  The  passage  of  the  shock  heats 
the  material  through  adiabatic  compression;  if  the  compression  is  great  enough,  then  the 
local  temperature  will  rise  to  a  threshold  required  to  initiate  the  chemical  reactions,  and  a 
reaction  front  will  form  behind  the  shock.  In  practice,  the  equations  of  state  for  energetic 
materials  are  quite  complex,  and  feature  many  empirical  parameters,  whose  adjustment 
allows  the  theory  to  agree  quite  well  with  preexisting  experimental  data. 

However,  for  understanding  the  detailed  microscopic  nature  of  detonations  in  a  solid, 
the  traditional  theories  may  not  be  appropriate.  At  the  scale  of  intermolecular  distances, 
the  usefulness  of  thermodynamic  variables  is  problematic,  and  at  timescales  shorter  than 
nanoseconds,  relations  built  upon  the  assumption  of  thermodynamic  equilibrium  require 
special  scrutiny  (the  rapid  approach  to  thermal  equilibrium  under  some  circumstances  is 
discussed  further  in  section  3.)  We  should  point  out  that  a  description  based  on  classical 
mechanics,  such  as  that  employed  here,  must  also  be  insufficient,  but  can  form  the  basis 
for  more  accurate  quantum  or  semi-classical  descriptions. 

In  this  paper  we  examine  the  point  of  view,  championed  by  Walker,4-'  that  a  shock 
in  a  solid  initiates  detonation  through  the  mechanical  generation  of  scission  forces  on  the 
molecules  comprising  the  solid,  breaking  chemical  bonds,  creating  a  distribution  of  free 
radicals,  and  supplying  the  kinetic  energy  required  to  initiate  reaction.  In  the  light  of  this 
picture  of  the  detonation  process,  we  investigate  the  interaction  of  the  shock  front  with 
the  lattice  structure,  using  numerical  molecular  dynamics.  We  are  particularly  interested 
in  the  effects  of  certain  types  of  lattice  defects  on  the  shock  initiation  mechanism.  We 
include  no  chemistry  in  our  model;  the  “molecules”  are  point  particles  interacting  through 
Lennard-Jones  potentials.  Hence  our  current  connection  to  the  detonation  problem  is  the 
exploration  of  how  the  conditions  leading  to  bond  scission  and  subsequent  recombination 
are  established.  The  following  stage,  of  reaction  and  shock  acceleration,  will  be  treated  in 
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a  subsequent  paper,  which  will  report  on  simulations  employing  a  chemistry  model. 

There  have  been  reported  a  number  of  interesting  molecular  dynamics  calculations  that 
bear,  in  one  way  or  another,  on  the  problem  of  the  shock  initiation  of  detonations  in  solids. 
Karo  et  al.s  simulated  a  small  lattice  of  approximately  one  hundred  molecules  arranged  in  a 
two-dimensional  square-symmetric  lattice,  interacting  through  either  “endothermic”  Morse 
potentials  or  “exothermic”  polynomial  potentials.  A  shock  was  launched  by  simulating  a 
flying  plate  consisting  of  a  smaller  lattice  of  identical  construction  to  the  main  lattice.  (This 
is  quite  similar  to  our  method  of  shock  launching  in  this  paper.)  The  main  conclusions  were 
that  the  shock  was  quite  narrow  on  the  atomic  scale,  that  the  temperature  of  the  lattice  was 
unimportant,  and  that  the  free  surface  at  the  end  of  the  finite  lattice  was  a  crucial  feature, 
as  the  shock’s  interaction  with  this  surface  resulted  in  a  spalling  off  of  a  large  chunk  of  the 
crystal.  The  behavior  subsequent  to  the  spalling  depended  on  which  potential  was  in  use. 
With  the  endothermic  potential,  there  was  little  interesting  activity  after  the  separation 
of  the  chunk,  but  with  the  exothermic  potential,  a  violent  reorganization  of  the  lattice 
occurred. 

While  there  is  little  reason  to  question  the  reliability  of  the  general  character  of  the 
results  reported  in  these  seminal  papers,  especially  the  importance  of  free  surfaces,  there 
are  several  numerical  issues  that  should  be  discussed.  Unfortunately,  issues  of  convergence 
and  accuracy  are  not  treated  in  these  references:  conserved  quantities  are  not  mentioned, 
nor  are  the  methods  of  integration.  Instead  of  physically  faithful  additive  potentials, 
the  authors  employ  two  different  first-  and  second-neighbor  potentials  chosen  to  make 
the  initial  lattice  stable.  In  addition,  as  the  lattice  undergoes  its  natural  distortion  and 
molecules  acquire  different  sets  of  close  neighbors,  they  are  not  allowed  to  interact  with 
these  new  neighbors.  Instead,  the  original  table  of  bonds  is  used  throughout  the  calculation. 
As  the  authors  point  out,  this  leads  at  times  to  molecules  passing  through  each  other 
without  interaction.  (Indeed,  this  is  a  statistical  possibility,  albeit  of  low  likelihood,9  with 
the  method  we  employ,  as  discussed  below.)  One  of  our  goals  in  the  work  reported  here 
is  to  discover  which  of  these  results  survive  the  application  of  more  modem  numerical 
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methods  applied  to  somewhat  more  realistic  simulations. 

Later,10  using  the  same  numerical  method,  these  authors  studied  a  shock  in  a  lattice 
interrupted  by  a  large  gap  consisting  of  a  region  with  no  molecules.  As  before,  they 
observed  spall  from  the  free  surface,  with  this  time  the  spalled  molecules  reaching  the 
second  free  surface  beyond  the  gap  and  launching  a  second  shock.  This  was  an  attempt 
to  simulate  efficiently  the  behavior  of  a  true  void  in  a  crystal,  which  would  be  surrounded 
by  the  lattice  on  all  sides  and  make  up  a  small  fraction  of  the  lattice.  As  we  discuss  in  the 
section  3,  our  simulations  of  small  true  voids  included  in  a  large  lattice  lead  to  somewhat 
different  conclusions  from  simulations  employing  a  gap. 

In  more  recent,  related  work,  these  authors  and  others11-13  embed  a  computational 
surrogate  for  a  polyatomic  molecule  into  a  simple  host  lattice,  and  show  how  the  passage 
of  the  shock  pumps  energy  directly  into  some  of  the  vibrational  modes  of  the  molecule. 
This  pumping  could  lead  directly  to  bond  scission  in  a  real  material,  in  a  region  that  is 
definitely  not  ergodic,  and  where  no  thermodynamic  temperature  can  be  defined. 

Also  relevant  is  the  work  of  Tsai  and  Trevino14  on  a  diatomic  perfect  crystal  in  the  form 
of  a  long  filament  with  a  small  cross-section.  Exothermic  bond  breaking  was  simulated 
through  the  use  of  bound  and  free  compound  Morse  potentials,  and  the  shock  was  formed 
by  heating  the  first  six  crystal  planes.  The  initial  temperature  of  their  lattice  was  set  at  just 
below  the  threshold  of  spontaneous  dissociation;  the  heating  caused  dissociation  directly, 
the  reaction  spread  by  thermal  conduction  and  drove  a  shock  wave  into  the  filament.  They 
kept  track  of  the  stresses  and  examined  the  partition  of  kinetic  energy  in  the  induction 
zone,  concluding  that  thermal  equilibrium  was  not  reached  before  dissociation  took  place. 
There  is  no  discussion  of  numerical  accuracy  or  the  treatment  of  distant  neighbors,  although 
presumably  the  expediency  of  imposing  a  cutoff  distance  on  the  potentials,  as  adopted  in 
a  previous  paper,15  was  used  here  as  well. 

Several  authors  have  demonstrated  that  the  typical  quasi-steady  detonation  structure, 
described  above,  is  predicted  by  molecular  dynamics  models,  with  the  exothermic  reactive 
process  modeled  either  by  classical  potentials  acting  between  diatomic  bonds16,9  or  by 
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a  more  elaborate  procedure  that  captures  some  of  the  features  of  the  quantum  chemical 
processes  involved.17-21  All  of  these  calculations  involve  perfect  crystals,  and  although  some 
of  them  examine  the  effects  of  lattice  geometry,  shock  direction  relative  to  lattice  symmetry, 
etc.,  none  of  them  deal  with  the  action  of  crystallographic  defects.  A  recent  exception  is 
the  work  of  Maffre  and  Peyrard,22  which  deals  with  the  ability  of  a  developed  detonation 
structure  to  traverse  grain  boundaries,  voids,  and  other  localized  defects.  However  they 
do  not  deal  with  the  effect  of  defects  in  the  transition  from  a  pure  shock  to  a  detonation. 
Another  exception  is  the  recent  work  of  Tsai,23  which  will  be  discussed  below. 

2  METHOD 

Our  model  system  consists  of  an  array  of  point  particles,  each  with  a  specified  initial 
position  and  velocity  vector,  arranged  in  a  two-dimensional  space  with  either  periodic 
or  free  boundary  conditions  in  both  directions.  The  algorithms  were  chosen  for  their 
convenience  and  efficiency  on  the  computer  used,  the  Connection  Machine  CM-200,  which 
consists  of  a  large  number  (4096,  8192,  or  16384,  depending  on  configuration)  of  processors, 
each  with  its  own  set  of  locally  stored  data,  updated  in  parallel  according  to  instructions 
from  a  controlling,  “front-end”  computer.  Efficiency  consists  here  largely  of  maximizing  the 
proportion  of  code  that  executes  in  parallel,  which  entails  minimizing  the  manipulation  of 
front-end  (global)  data  and  maintaining  careful  control  over  the  pattern  of  communication 
among  processors. 

Forces  between  particles  were  represented  by  Lennard-Jones  potentials.  The  solution 
of  the  complete  N-body  problem  was  not  attempted;  rather,  the  short-range  nature  of  the 
Lennard-Jones  potential  was  exploited  to  limit  the  distance  over  which  interactions  had 
to  be  computed.  The  particle  data  was  distributed  onto  a  data  structure  known  as  the 
Monotonic  Lagrangian  Grid  (MLG).24  This  is  an  object  tracking  and  sorting  technique 
where  each  particle  is  associated  with  a  pair  of  integer  indices  i  and  j,  and  the  assignment 
is  ordered  such  that  the  coordinate  x  increases  monotonically  with  t,  and  y  increases 
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monotonically  with  j.  Monotonicity  was  enforced  after  each  update  in  particle  positions 
by  a  parallel  version  of  the  swapping  routine  described  in  reference  24.  The  advantage 
of  this  method  is  that,  at  each  timestep,  a  particle’s  neighbors  can  be  identified  quickly 
by  cycling  over  a  predetermined  region  of  i  —  j  space,  rather  than  by  searching  through 
coordinate  space  with  an  expensive  comparison  of  particle  separations. 

Each  physical  processor  on  the  Connection  Machine  can  be  divided,  in  software,  into  as 
many  virtual  procesors  as  memory  will  allow.  When  we  refer  to  a  “processor”  from  here  on, 
this  should  be  taken  to  mean  “virtual  processor”.  One  particle  was  associated  with  each 
processor,  which  stored  the  particle’s  position  and  velocity  coordinates,  as  well  as  any  other 
information  unique  to  the  particle.  We  also  stored  a  near  neighbors  template24  in  each 
processor;  this  is  an  array  that  stores  the  position  information  for  each  particle’s  neighbors 
in  grid-index  space.  Once  the  template  is  filled,  the  total  force  vector  acting  on  each 
particle  and  the  update  in  position  is  computed  for  all  particles  entirely  in  parallel,  with 
no  additional  interprocessor  communication  needed  for  the  rest  of  the  timestep.  After 
each  regeneration  of  the  MLG  by  the  swapping  routine,  the  template  is  rebuilt;  thus 
each  processor-particle  always  has  a  locally  stored,  updated  list  of  the  coordinates  of  its 
near  neighbors.  The  use  of  the  MLG  is  one  of  several  available  efficient  schemes  for 
replacing  a  time  consuming  search  for  the  spatially  proximate  neighbors  of  each  particle  at 
every  timestep.  However,  as  alluded  to  in  the  Introduction,  its  use  can  lead  occasionally 
(but  infrequently)  to  the  problem  of  missed  near  neighbors,24  particles  that  are  spatially 
close  but  not  included  in  the  neighbor  template.  When  these  near-misses  are  detected 
by  a  failure  of  energy  conservation  or  repeatability,  we  redo  the  calculation  with  a  larger 
template.  Thus  the  template  size  is  treated  like  the  computational  timestep,  both  of  which 
are  adjusted  to  keep  the  calculations  stable  and  accurate. 

The  procedure  for  filling  up  the  template  at  each  timestep  is  an  important  determinant 
of  the  overall  efficiency  of  the  algorithm.  Interprocessor  communication  between  neighbor¬ 
ing  processors  happens  to  be  far  more  efficient  than  between  more  distant  processors,  so 
we  fill  the  template  using  a  sequence  of  coordinated  data  moves  between  neighboring  pro- 
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cessors  only.  The  total  number  of  these  parallel  move  operations  is  equal  to  one  less  than 
the  total  number  of  elements  in  the  template  array;  therefore  the  total  time  required  to  fill 
the  template  depends  on  its  size,  not  on  the  number  of  particles  (assuming  that  the  ratio 
of  virtual  processors  to  physical  processors  remains  constant).  Part  of  the  communication 
sequence  is  diagrammed  schematically  in  figure  1. 

As  the  studies  described  in  this  paper  were  being  carried  out,  the  computer  hardware 
and  system  software  being  used  was  evolving.  Fortunately,  only  one  of  these  changes  needs 
to  be  discussed  here:  initially,  the  C  jnnection  Machines  that  were  used  for  these  calcula¬ 
tions  were  equipped  with  floating-point  accelerators  that  operated  on  32-bit  numbers  only; 
more  recently,  we  were  provided  with  64-bit  accelerators. 

For  reasonable  efficiency,  it  is  necessary  to  restrict  the  precision  of  floating  point  num¬ 
bers  to  the  word  size  of  the  accelerator.  Therefore,  some  of  the  results  reported  here  were 
performed  with  32  bit  arithmetic,  some  with  64  bit,  and  some  of  the  32  bit  calculations 
were  redone  with  64  bits  (“double  precision”)  when  that  b^  ame  available.  The  small  dif¬ 
ferences  we  observed  in  particle  data  between  the  two  precisions  was  not  significant  enough 
to  affect,  our  conclusions.  The  chief  advantage  to  the  greater  precision  will  be  for  future 
calculations,  which  we  will  be  able  to  carry  out  for  longer  times  before  the  accumulated 
errors  become  unacceptable. 

We  found  that  single  precision  arithmetic  limited  the  accuracy  with  which  energy  can 
be  conserved,  and  made  the  extra  accuracy  of  higher-order  time  integration  methods  su¬ 
perfluous,  and  their  use  needlessly  time  consuming.  Therefore  some  of  the  results  reported 
here  were  calculated  with  the  simple  explicit  Euler  method.  We  were  able  to  acheive  sta¬ 
bility  with  this  method  using  a  reasonable  timestep.  The  consequence  of  the  limitation 
of  precision  is  that  energy  is  conserved  only  to  approximately  one  percent  in  the  32  bit 
calculations.  With  double  precision  arithmetic,  it  became  advantageous  to  use  a  higher- 
order  integrator.  We  have  used  both  the  leapfrog  method  (described  in  several  places,  for 
instance25)  and  the  third  order  Verlet  formula26  combined  with  a  velocity  corrector  de¬ 
scribed  by  Beeman.27  The  advantage  of  this  method  (described  at  the  bottom  of  page  138 
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of  reference  27)  is  that  the  velocities  do  not  have  to  be  permanently  stored,  and  need  only 
be  calculated  when  they  are  required  for  a  data  dump  or  a  check  of  the  kinetic  energy. 
Our  algorithm  with  these  integrators  can  conserve  energy  to  one  part  in  104,  in  actual 
calculations. 

We  have  performed  calculations  with  three  different  system  sizes,  “small”  grids  of  212  or 
213  molecules  distributed  among  an  equal  number  of  physical  processors,  and  a  “large”  grid 
of  215  molecules  distributed  among  214  physical  processors,  which  is  the  largest  processor 
set  available  to  us.  (Our  original  machine,  with  two  banks  of  212  processors,  was  replaced 
partway  through  thic  project  with  one  of  two  banks  of  213  processors.)  In  both  cases  each 
molecule  was  assigned  its  own  virtual  processor.  All  else  being  equal,  the  large  grid  should 
take  approximately  twice  as  long  to  simulate  as  the  small  grids,  because  it  is  the  ratio  of 
virtual  to  physical  processors  that  is  significant.  On  the  small  grids,  with  five  neighbors 
in  the  template  in  each  direction,  the  calculation  took  0.67  seconds  per  timestep.  On  the 
large  grid,  a  five-neighbor  calculation  took  1.0  seconds  per  timestep  and  a  seven-neighbor 
calculation  took  2.0  seconds  per  timestep.  We  have  verified  several  of  our  runs  with 
equivalent  calculations  on  NRL’s  Cray  X-MP,  using  a  well-vectorized  code.  With  4096 
molecules,  a  five-neighbor  template  takes  2-2.4  seconds  per  timestep;  with  larger  system 
sizes,  the  disparity  between  the  two  architectures  can  be  expected  to  grow  dramatically.28 

3  RESULTS 

It  is  difficult  to  compare  directly  the  various  molecular  dynamics  simulations  that  can  be 
found  in  the  literature,  due  to  the  varying  sets  of  units  employed  by  different  authors, 
along  with  differences  in  potential  parameters,  lattice  spacing,  initial  conditions,  etc.,  each 
of  which  can  have  a  nonobvious  effect  on  the  evolution  of  the  system.  Although  the 
large  number  of  parameters  characterizing  a  molecular  system  makes  a  simple  resort  to 
dimensionless  numbers  impossible,  it  may  be  useful  to  introduce  one  such  number,  which 
can  be  used  to  organize  the  results  concerning  the  type  of  molecular  dynamics  configuration 
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used  in  this  paper. 

We  have  an  array  of  molecules  with  forces  between  them  derived  from  a  Leonard- J ones 
potential, 

<t>  =  4c  [(<r/r)12  -  (ff/r)6]  , 

so  the  force,  F,  is  given  by 


dd> 


<Pr 
mdi *' 


by  Newton’s  equation  of  motion.  In  the  above  two  equations  r  is  the  distance  from  the 
molecule,  m  is  the  molecular  mass,  t  is  time,  an  and  6  are  the  two  parameters  char¬ 
acterizing  the  potential.  In  order  to  make  the  .ation  of  motion  dimensionless,  we  must 
introduce  scales  appropriate  to  a  system  excited  by  a  flying  plate: 


d  =  interparticle  spacing 

fi  —  particle  mass 
After  making  the  substitutions 

r  r'd  e  —  e'E 


Vp  =  plate  velocity  r  =  d/Vp 
E  =  fid2 jr7  =  energy  unit. 

t  —*  t'r  m  — »  m'fi  <j  — ♦  a'd 


and  removing  the  primes  from  the  now  dimensionless  variables,  the  equation  of  motion 
becomes 


£-'4  [— 12<r,J/r,3  +  6^/rr]  =m^J, 


dt2 


where 


emerges  as  a  new  dimensionless  number,  apparently  relating  the  kinetic  energy  of  the 
impacting  molecules  to  the  binding  energy  of  the  intermolecular  potential. 

In  all  of  the  calculations  presented  here,  e  =  0.0223,  a  =  0.891,  the  plate  velocity  = 
2.0,  the  intermolecular  separation  at  equilibrium  =  1.0,  and  the  particle  mass  =  4.0,  giving 


10 


Phillips  et  al .:  Molecular  Dynamics  of  Shocks.... 


C  =  727.3.  It  may  help  put  things  into  context  for  some  readers  if  we  scale  our  variables 
to  the  values  appropriate  to  solid  argon,  a  substance  commonly  discussed.  In  this  case  the 
plate  velocity  scales  to  4258  meters/second,  one  time  unit  is  0.88  x  10-13  seconds,  and  one 
distance  unit  is  3.8  A.  From  here  on  we  shall  use  the  dimensionless  units. 

Figure  2  shows  a  sequence  of  molecular  configurations  as  a  plate-launched  shock  en¬ 
counters  a  pair  of  voids,  producing  a  great  deal  of  disorder.  The  same  initial  conditions 
in  a  system  without  the  voids  lead  to  the  shock  traversing  the  crystal  intact  and  leaving 
it  undisturbed.  We  see  here  that  the  voids  disturb  the  coherency  of  the  shock,  trans¬ 
forming  its  organized,  x-directed  (horizontal)  motions  into  a  thermalized,  two-dimensional 
(x  —  y)  motion;  beyond  the  defects,  the  shock  continues  as  several  disconnected  pieces.  The 
progress  of  the  shock  can  perhaps  be  more  clearly  visualized  in  figure  3,  showing  profiles 
of  the  horizontal  velocity  along  different  lines  through  the  crystal.  That  these  disordered 
regions  are  actually  “thermalized”  can  be  seen  by  looking  at  the  speed  distribution  of  the 
molecules  in  the  crystal.  Figure  4  is  a  series  of  speed  histograms  corresponding  to  the  cal¬ 
culation  illustrated  in  figure  2;  superimposed  on  each  histogram  plot  is  the  two-dimensional 
Maxwell-Boltzman  distribution  with  a  temperature  derived  from  the  average  kinetic  en¬ 
ergy  of  all  the  particles.  At  early  times  we  can  see  the  speed  distribution  dominated  by  two 
values,  the  near  zero  speed  of  the  particles  locked  into  the  near  zero  temperature  lattice 
configuration,  and  the  highest  speed,  equal  to  the  plate  velocity.  As  time  advances  we  can 
see  the  velocity  distribution  filling  in,  but  remaining  non- Maxwellian.  In  the  last  frame 
of  the  figure  we  have  plotted  the  distribution  of  a  subset  of  the  system  at  t  =  15  (the 
configuration  at  this  time  is  shown  in  Figure  15),  the  subset  consisting  of  the  particles  in 
the  chaotic  region  delimited  by  x  =  15  and  x  =  42.  The  superimposed  Maxwell-Boltzman 
curve  corresponds  to  the  kinetic  temperature  of  the  subset  of  particles  under  consideration; 
clearly  this  region  of  particles  has  rapidly  attained  a  thermal  equilibrium,  indicating  that 
the  collisionality  in  the  chaotic  region  is  high. 

A  similar  calculation  can  be  seen  in  figure  5,  with  the  square  lattice  replaced  by  an 
hexagonal  lattice.  In  this  case,  also,  the  small  void  has  been  replaced  by  a  single  va¬ 
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cancy.  We  can  see  here  that  even  in  a  lattice  geometry  with  a  higher  binding  energy,  and 
with,  unlike  the  square  lattice,  nonlinear  stability  against  moderate  temperatures,  a  small 
imperfection  creates  a  large  disruption. 

The  creation  of  disordered  hot  spots  from  shock- void  interaction  has  also  been  seen  in 
recent  simulations  by  Tsai.23  These  disordered  regions  also  have  the  property  of  being  at 
a  slightly  higher  density  than  the  equilibrium  lattice,  as  can  be  seen  in  the  density  plot  of 
Figure  6. 

The  interaction  of  a  shock  front  with  a  mass  defect  in  square  and  hexagonal  lattices 
can  be  seen  in  figures  7  and  8,  both  showing  the  effect  of  the  inclusion  of  a  molecule  with 
a  %50  mass  excess.  The  effects  are  similar  to  but  weaker  than  those  caused  by  voids. 

We  have  performed  a  series  of  experiments  on  systems  identical  to  the  ones  shown 
in  the  figures,  but  at  various  temperatures  and  intermolecular  potential  well  depths.  At 
very  small  temperatures  the  perfect  square  lattice  becomes  unstable  when  struck  by  the 
plate,  and  begins  to  dissociate.  This  is  true  even  if  the  potential  well  depth  is  increased 
by  a  factor  of  45.  The  hexagonal  lattice,  however,  is  stable  to  impact  by  the  flying  plate 
over  a  wide  range  of  temperatures,  although  at  higher  temperatures  the  shock  front  is 
less  sharply  defined.  The  effects  of  the  various  types  of  hexagonal  lattice  defects,  and 
their  relative  importance,  is  not  affected  by  temperature,  but  is  strongly  dependent  on  the 
intermolecular  potential  well  depth.  For  a  given  shock  strength  and  defect  type,  it  appears 
that  it  is  possible  to  increase  the  well  depth  to  a  point  at  which  the  lattice  remains  stable. 

It  is  useful  to  quantify  the  concept  of  the  strength  of  the  effect  of  different  types  of 
lattice  defects,  or  their  efficacy  in  disrupting  the  crystal  when  interacting  with  a  shock. 
For  this  purpose  we  have  defined  two  “disruption  factors”,  to  be  used  with  the  two  lattice 
symmetries.  The  disruption  factors  should  be  defined  in  such  a  way  that  the  elastic  defor¬ 
mation  of  the  lattice  due  to  the  shock  itself  does  not  make  a  contribution.  Therefore  the 
disruption  factors  will  be  zero  for  the  case  of  a  shock  traversing  a  defect-free  crystal  that 
remains  intact. 

We  can  define  disruption  factors  that  satisfy  these  criteria  as  follows.  For  the  square 
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lattice  the  angular  positions  of  the  four  nearest  neighbors  are  found.  The  disruption  factor, 
X,  for  each  particle  is  defined  as 

4 

X  =  X)  I  sin  0,  cos  0,1, 

»=i 

where  the  summation  is  over  near  neighbors.  For  a  compression  or  rarefaction  of  the 
lattice  due  to  the  propagation  of  a  shock  along  one  of  the  principal  crystal  axes,  the 
angular  positions  of  the  near  neighbors  should  remain  0,  ±tt/2,  and  7r  resulting  in  a  x  =  0. 

The  hexagonal  lattice  disruption  factor  must  be  defined  in  a  slightly  different  manner. 
The  angular  positions  of  the  six  near  neighbors  of  each  particle  change  as  the  shock  prop¬ 
agates  through  the  lattice.  In  this  case  the  angular  positions  of  the  six  near  neighbors  are 
first  found.  The  neighbors  in  the  lower  half-plane  are  reflected  into  the  upper  half-plane. 
The  near  neighbors  are  then  sorted  into  order  of  increasing  angle.  The  disruption  factor 
for  the  hexagonal  lattice  is  defined  as 

X  =  sin  0!  +  sin(03  —  02)  +  sin(05  —  0«)  +  sin  06. 

In  the  case  of  a  shock  passing  through  a  hexagonal  lattice  causing  only  a  compression  or 
rarefaction  along  the  shock  direction,  0i  and  06  should  remain  0  and  it,  respectively.  The 
angles  03  and  02  should  be  equal  and  05  and  04  should  be  equal.  As  in  the  square  lattice, 
a  zero  disruption  is  calculated. 

The  total  disruption  of  the  lattice  is  calculated  by  summing  the  per  particle  disruption 
factor  over  all  particles.  The  values  of  the  total  disruption  correlate  well  with  the  apparent 
disruption  found  by  visual  inspection  of  the  particle  positions. 

Figures  9  and  10  display  x  as  a  function  of  time  for  an  assortment  of  crystal  defects, 
for  the  square  and  hexagonal  lattices  respectively.  The  greater  importance  of  voids  over 
mass  defects  is  obvious,  as  is  the  unstable  growth  of  the  disordered  region  after  the  passage 
of  the  shock.  Figure  9  also  shows  that  the  fastest  growing  disruption  factor  is  found  for 
an  interstitial  inclusion  in  a  square  lattice.  The  spatial  configuration  for  such  a  system  is 
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shown  in  figure  11.  The  closer  packed  hexagonal  lattice  does  not  provide  an  equilibrium 
position  for  an  interstitial  inclusion. 

Another  useful  diagnostic  summary  of  the  state  of  the  system  is  a  histogram  distribution 
of  the  near-neighbor  angular  positions,  such  as  that  shown  in  figure  12.  This  figure  shows 
the  evolution  of  the  angular  distribution  for  the  case  of  a  vacancy  in  a  square  lattice.  The 
additional  peaks  which  develop  in  the  histogram  distributions  indicate  that  the  lattice 
is  beginning  to  locally  undergo  a  transition  from  a  square  to  a  hexagonal  state.  For  a 
hexagonal  lattice  formed  by  the  displacement  of  every  other  column  of  atoms  perpendicular 
to  the  shock  direction  together  with  a  compression  along  the  shock  direction,  the  near 
neighbor  angular  positions  should  be  ±tt/6,  ±it/2  and  ±5ir/6.  The  peaks  are  not  found 
exactly  at  these  angles  since  the  lattice  has  a  square  structure  a  short  distance  from  the 
region  of  disturbance,  resulting  in  a  distorted  hexagonal  state. 

Figure  13  shows  a  calculation  similar  to  that  shown  in  figure  2,  but  in  a  system  of 
32K  molecules,  with  a  distribution  of  voids  whose  positions,  shapes,  and  sizes  (ranging 
from  single-site  vacancies  to  voids  of  nine  molecules)  are  determined  by  a  random  number 
generator.  This  is  a  first  attempt  at  simulating  a  system  approximating  a  part  of  a  real 
crystal,  with  its  natural  distribution  of  defects.  The  basic  process  shown  in  figure  2  is 
repeated  here:  each  void,  including  even  the  single-site  vacancies,  becomes  the  seed  of  a 
rapidly  growing  region  of  thermalized  disorder,  which  blocks  the  progress  of  the  shock. 
The  shock  becomes  more  tenuous  as  its  organized  i-momentum  is  equilibrated,  and,  in 
the  absence  of  energy  available  from  chemical  reactions,  eventually  dies  out. 

It  is  of  some  relevance  to  the  mechanical  bond-scission  view  of  detonation  initiation4 
to  determine  the  character  of  the  forces  acting  upon  individual  molecules.  While  the  effect 
of  intermolecular  potentials  on  the  breaking  of  bonds  is  complicated  and  difficult  to  treat 
faithfully  with  a  classical  description,  we  believe  we  can  make  some  relevant  observations. 
The  quantities  of  interest  are  the  components  of  the  stress  tensor  at  each  molecule;  however, 
since  our  “molecules”  are  point  particles,  with  no  preferred  directions  for  bond  breaking, 
a  dilatory  force  along  a  certain  axis  combined  with  a  compressional  force  along  a  different 
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axis  can  not  be  interpreted  as  an  influence  that  would  necessarily  lead  to  the  breaking  of 
a  bond.  Our  compromise  is  to  look  at  the  divergence  of  the  force  vector,  V  •  F,  a  quantity 
that  when  positive  indicates  a  net  extensional  force  acting  on  the  molecule.  Figure  14  is  a 
particle  position  plot  corresponding  to  a  frame  of  figure  2,  but  indicating  which  molecules 
are  experiencing  a  positive  force  divergence  and  which  are  experiencing  a  negative  VF.  We 
see  from  this  figure  that  the  molecules  participating  in  the  advance  of  the  shock  front  are 
under  a  large  compressional  force,  while  some  of  the  molecules  in  the  disordered  regions 
behind  the  shock  have  V  •  F  >  0,  and  so  are  being  “pulled  apart”  by  their  neighbors. 
However,  the  magnitude  of  V  •  F  for  these  molecules  is  much  smaller  than  that  for  the 
compressed  molecules  in  the  shock  front. 

The  spherically  symmetric  potential  leads  to  two  possible  periodic  equilibrium  lattice 
symmetries:  square  and  hexagonal.  As  can  be  seen  in  Figure  15,  the  agitation  created  by 
the  void  collapse  has  provided  a  section  of  the  lattice  the  opportunity  to  relax  to  the  other 
periodic  state  available  to  it:  we  see  here  and  in  various  other  rims  the  emergence  of  a 
localized  hexagonal  region.  We  use  the  term  “relax”  above  because  the  hexagonal  state  is 
actually  of  slightly  smaller  potential  energy.  Accompanying  this  relaxation,  therefore,  must 
be  the  generation  of  some  extra  kinetic  energy.  A  somewhat  similar  shock-induced  phase 
transition  was  shown  in  reference  10  where  a  potential  with  two  minima  acted  between 
molecules. 

The  assertion  has  been  made10  that  in  order  to  understand  the  shock-void  interaction 
it  is  not  necessary  to  place  the  void  interior  to  the  crystal,  but,  in  the  interests  of  com¬ 
putational  efficiency  (which  in  this  context  means  reducing  the  number  of  molecules  that 
need  to  be  tracked)  it  is  sufficient  to  examine  the  interaction  of  the  shock  with  a  gap  in 
the  lattice.  In  figure  16  we  show  the  results  of  a  calculation  in  every  way  identical  with 
that  shown  in  figure  2,  except  that  the  two  voids  in  figure  2  have  been  replaced  with  a  gap 
extending  across  the  entire  crystal;  the  gap  begins  at  the  same  x  position  as  the  leading 
edges  of  the  voids,  and  has  the  same  width.  We  can  see  from  Figure  16  that  the  shock  suc¬ 
ceeds  in  traversing  the  gap  while  maintaining  its  coherency,  and  proceeds  down  the  lattice, 
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leaving  it  undisturbed.  Thus  we  see  an  essential  physical  difference  between  a  gap  and  an 
included  void:  the  presence  of  neighbor  molecules  transverse  to  the  direction  of  progress 
of  the  shock  causes  energy  to  be  removed  from  the  shock  and  distributed  among  greater 
degrees  of  freedom,  leading  to  a  thermalized  region.  A  related  observation  has  recently 
been  made  by  Maffre  and  Peyrard,22  who  demonstrated  the  ability  of  a  shock-reaction 
front  structure  to  traverse  a  void  in  a  lattice. 

In  addition  to  various  systems  of  voids,  we  have  simulated  the  encounter  of  a  shock 
with  a  type  of  extended  defect  sometimes  called  a  slip  line;  the  evolution  of  the  system  is 
shown  in  figure  17,  which  contains  a  sequence  of  plots  of  the  particle  positions.  The  slip 
line  can  be  seen  at  x  =  24.5;  there  are  also  small  voids  at  the  top  and  bottom  boundaries 
extending  from  the  slip  line,  created  by  the  periodic  boundary  conditions.  Once  again  the 
voids  are  seeds  of  a  single  growing  chaotic  region  (joined  at  the  periodic  boundary),  while 
it  is  clear,  by  inspecting  the  central  portion  of  the  plots,  that  the  shock  traverses  the  slip 
line  completely  intact. 


4  CONCLUSIONS 

When  a  shock,  sufficiently  weak  that  it  is  able  to  traverse  a  perfect  crystal  without  per¬ 
manently  disturbing  the  configuration  of  its  lattice,  encounters  a  void  in  the  lattice,  the 
void  becomes  the  site  of  a  rapidly  growing  thermalized,  hot,  fluid-like  phase  character¬ 
ized  by  a  high  density  and  a  high  degree  of  collisionality.  These  characteristics  should  be 
conducive  to  the  onset  of  chemical  reactions,  and,  we  suspect,  it  is  in  these  regions  that 
the  reactions  leading  to  the  development  of  a  shock-detonation  stucture  begin.  It  seems 
probable,  therefore,  that  the  void  distribution  in  the  lattice  is  an  important  factor  con¬ 
trolling  the  sensitivity  to  shock-initiation  and  the  character  of  the  subsequent  detonation 
front  development.  A  perfect  crystal  should  be  relatively  insensitive.  It  is  possible  that  in 
three  dimensions,  other  types  of  defects  will  be  seen  to  be  equally  important,  but  that  will 
be  treated  in  a  subsequent  paper. 


16 


Phillips  ct  al .:  Molecular  Dynamics  of  Shocks.... 


As  discussed  above,  in  order  to  concentrate  on  the  narrowly  defined  problem  of  shocks 
and  defects  in  molecular  crystals,  with  as  few  complications  as  possible,  we  have  deferred 
including  a  model  of  the  chemical  bond  and  simulated  the  behavior  of  moderately  large 
systems  of  indivisible  molecules  with  spherically  symmetrical  potentials.  This  approach 
has  the  advantage  that  our  results  axe  independent  of  the  details  of  any  particular  model 
of  detonation  chemistry  or  intramolecular  behavior,  several  of  which  are  referenced  in  the 
Introduction. 

Although  we  have  been  able  to  resolve  several  issues  concerning  the  importance  of 
defects,  there  are  some  remaining  ambiguities  that  will  not  be  resolved  until  we  include 
polyatomic  molecules  in  the  simulations.  These  have  their  counterpart  in  the  uncertainties 
plagueing  our  knowledge  of  the  bond  scission  process  in  a  shocked  crystal  lattice.  If  the 
polyatomic  bonds  are  broken  in  the  shocked  region,  due  to  direct  energy  transfer  from  the 
shock  to  vibrational  modes  of  the  molecules,  then  the  scission  forces  discussed  above  have 
little  relevance,  because  they  occur  in  the  disordered  regions  behind  the  shock.  In  this 
case  the  importance  of  the  defects  is  in  the  thermalization  and  mixing,  which  will  provide 
enhanced  opportunities  for  the  free  radicals,  created  in  the  shocked  region,  to  recombine. 
Of  particular  relevance  here  is  the  observation  that  thermal  equilibrium  is  established 
in  the  disordered  region  close  behind  the  shock,  on  a  very  fast  timescale,  implying  that 
equilibrium  equations  of  state  may  be  relevant  after  all  to  the  detonation  process.  If  the 
chemical  bonds  are  not  broken  directly  by  the  passing  shock,  then  the  scission  forces  may 
be  responsible  for  bond  stretching  and  breaking  in  the  disordered  region,  where  conditions 
prevail  that  will  enhance  the  subsequent  reactivity.  The  truth  is  that  neither  computational 
nor  theoretical  work  to  date  is  sufficient  to  resolve  these  uncertainties.  We  hope  that  these 
questions  can  be  addressed  in  the  next  generation  of  simulations,  involving  both  defects  and 
chemical  reactions  in  polyatomic  crystal  lattices  large  enough  to  capture  their  interaction. 
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Figure  1:  The  two  initial  steps  in  a  sequence  of  parallel  move  operations; 
the  simplified  case  of  a  two-neighbor  template  is  shown  as  an  illustrative 
example,  and  the  templates  for  six  of  the  particles  are  included.  Each 
processor’s  template  is  seeded  at  its  center  with  its  own  particle’s  data;  the 
arrows  represent  the  copying  of  particle  data  from  a  template  entry  to  the 
appropriate  template  entry  in  a  neighboring  processor. 
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Figure  2:  Particle  positions  for  4k  particles  in  two  dimensions.  A  “flying 
plate”  impinges  from  the  left,  and  has  just  made  contact  at  t  =  3;  the 
resulting  shock  interacts  with  two  rectangular  voids. 
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Figure  2:  Particle  positions  for  4k  particles  in  two  dimensions.  A  “flying 
plate”  impinges  from  the  left,  and  has  just  made  contact  at  t  —  3;  the 
resulting  shock  interacts  with  two  rectangular  voids. 
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Figure  4:  Velocity  distribution  histograms  with  superimposed  equilibrium 
distributions,  for  the  same  calculation  illustrated  in  figure  2. 


Figure  5:  A  void  in  a  hexagonal  lattice  interacting  with  a  shock  launched 
by  a  flying  plate. 
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Figure  6:  Scaled,  smoothed  density  plot  corresponding  to  the  same  calcu¬ 
lation  illustrated  in  figure  2,  at  t  —  13.5.  Increasing  brightness  indicates 
increasing  density.  The  three  bright  patches  near  the  right  hand  side  of  the 
figure  correspond  to  the  compression  at  the  location  of  the  void- interrupted 
shock. 
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Figure  7:  A  mass  defect  in  a  square  lattice  interacting  with  a  shock  launched 
by  a  flying  plate.  The  black  square  is  a  molecule  with  a  mass  excess  of  %50. 
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Figure  8:  A  mass  defect  in  a  hexagonal  lattice  interacting  with  a  shock 
launched  by  a  flying  plate.  The  black  square  is  a  molecule  with  a  mass 
excess  of  %50. 
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POSITIONS  AT  TIME  =  7.500E+00. 
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Figure  11:  An  interstitial  inclusion  in  a  square  lattice. 
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Figure  12:  The  evolution  of  the  distribution  of  intennolelcular  bond  angles 
in  a  square  lattice  as  the  shock  front  passes  through  a  lattice  vacancy. 
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PARTICLE  POSITIONS  AT  TIME  -0 


PARTICLE  POSITIONS  AT  TIME  =  1 . 101E-01 
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Figure  13:  Particle  positions  in  a  32k  particle  system,  with  a  random  dis¬ 
tribution  of  voids.  Again  a  flying  plate  has  launched  a  shock  from  the 
left.  Only  the  left  half  (approximately  16k  particles)  of  the  total  system  is 
shown. 


Figure  14:  The  quantity  V  •  F  for  calculation  illustrated  in  figure  2,  at 
t  =  10.5.  “+”  symbols  indicate  the  positions  of  particles  with  a  V  •  F  >  0, 
and  “M”  indicates  the  positions  of  particles  with  a  V  •  F  <  0. 
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PARTICLE  POSITIONS  AT  TIME  =  1.500E+01 


Figure  15:  A  frame  from  a  later  time,  t  =  15,  in  the  calculation  illustrated 
in  figure  2,  showing  the  emergence  of  an  hexagonal  phase  embedded  in  the 
predominant  square-symmetric  crystal.  A  box  has  been  drawn  around  the 
hexagonal  region. 
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Figure  16:  A  sequence  of  particle  position  plots  similar  to  those  in  figure  2, 
but  with  the  two  voids  replaced  by  a  continuous  gap. 
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Figure  16:  A  sequence  of  particle  position  plots  similar  to  those  in  figure  2, 
but  with  the  two  voids  replaced  by  a  continuous  gap. 
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Figure  17:  Particle  position  plot  of  a  shock  progressing  through  a  linear 
defect  and  two  small  voids  at  x  =  24.5 
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